import cartopy.crs as ccrs
import intake
import uxarray as ux
cat_url = "https://digital-earths-global-hackathon.github.io/catalog/catalog.yaml"
cat = intake.open_catalog(cat_url).NCAR
cat
NCAR:
  args:
    path: https://digital-earths-global-hackathon.github.io/catalog/NCAR/catalog.yaml
  description: catalog as visible from NCAR
  driver: intake.catalog.local.YAMLFileCatalog
  metadata:
    catalog_dir: https://digital-earths-global-hackathon.github.io/catalog
# cat.walk(depth=-1)
list(cat)
['CERES_EBAF',
 'ERA5',
 'IR_IMERG',
 'JRA3Q',
 'MERRA2',
 'casesm2_10km_nocumulus',
 'icon_d3hp003',
 'icon_d3hp003aug',
 'icon_d3hp003feb',
 'icon_ngc4008',
 'ifs_tco3999-ng5_deepoff',
 'ifs_tco3999-ng5_rcbmf',
 'ifs_tco3999-ng5_rcbmf_cf',
 'mpas_dyamond1',
 'mpas_dyamond2',
 'mpas_dyamond3',
 'nicam_gl11',
 'scream-dkrz',
 'scream2D_hrly',
 'scream_lnd',
 'scream_ne120',
 'tracking-d3hp003',
 'um_Africa_km4p4_RAL3P3_n1280_GAL9_nest',
 'um_CTC_km4p4_RAL3P3_n1280_GAL9_nest',
 'um_SAmer_km4p4_RAL3P3_n1280_GAL9_nest',
 'um_SEA_km4p4_RAL3P3_n1280_GAL9_nest',
 'um_glm_n1280_CoMA9_TBv1p2',
 'um_glm_n1280_GAL9',
 'um_glm_n2560_RAL3p3',
 'wrf_conus',
 'wrf_samerica']
ds_mpas2 = cat.mpas_dyamond2(zoom=3).to_dask()
ds_mpas2
/glade/u/apps/opt/conda/envs/2025-digital-earths-global-hackathon/lib/python3.12/site-packages/intake_xarray/base.py:21: FutureWarning: The return type of `Dataset.dims` will be changed to return a set of dimension names in future, in order to be more consistent with `DataArray.dims`. To access a mapping from dimension names to lengths, please use `Dataset.sizes`.
  'dims': dict(self._ds.dims),
<xarray.Dataset> Size: 357MB
Dimensions:                 (time: 329, cell: 768, nVertLevels: 55,
                             nSoilLevels: 4, nVertLevelsP1: 56)
Coordinates:
  * cell                    (cell) float64 6kB nan 1.0 2.0 ... 765.0 766.0 767.0
  * nSoilLevels             (nSoilLevels) float64 32B nan 1.0 2.0 3.0
  * nVertLevels             (nVertLevels) float64 440B nan 1.0 2.0 ... 53.0 54.0
  * nVertLevelsP1           (nVertLevelsP1) float64 448B nan 1.0 ... 54.0 55.0
  * time                    (time) datetime64[ns] 3kB 2020-01-20 ... 2020-03-01
Data variables: (12/19)
    air_temperature         (time, cell, nVertLevels) float32 56MB dask.array<chunksize=(1, 768, 55), meta=np.ndarray>
    h_oml                   (time, cell) float32 1MB dask.array<chunksize=(1, 768), meta=np.ndarray>
    hu_oml                  (time, cell) float32 1MB dask.array<chunksize=(1, 768), meta=np.ndarray>
    hv_oml                  (time, cell) float32 1MB dask.array<chunksize=(1, 768), meta=np.ndarray>
    pressure                (time, cell, nVertLevels) float32 56MB dask.array<chunksize=(1, 768, 55), meta=np.ndarray>
    relhum                  (time, cell, nVertLevels) float32 56MB dask.array<chunksize=(1, 768, 55), meta=np.ndarray>
    ...                      ...
    tslb                    (time, cell, nSoilLevels) float32 4MB dask.array<chunksize=(1, 768, 4), meta=np.ndarray>
    uReconstructMeridional  (time, cell, nVertLevels) float32 56MB dask.array<chunksize=(1, 768, 55), meta=np.ndarray>
    uReconstructZonal       (time, cell, nVertLevels) float32 56MB dask.array<chunksize=(1, 768, 55), meta=np.ndarray>
    vegfra                  (time, cell) float32 1MB dask.array<chunksize=(1, 768), meta=np.ndarray>
    w                       (time, cell, nVertLevelsP1) float32 57MB dask.array<chunksize=(1, 768, 56), meta=np.ndarray>
    xice                    (time, cell) float32 1MB dask.array<chunksize=(1, 768), meta=np.ndarray>

Let us plot snow

# Open a HEALPix UX array
uxds = ux.UxDataset.from_healpix(ds_mpas2)
uxds
<xarray.UxDataset> Size: 357MB
Dimensions:                 (time: 329, n_face: 768, nVertLevels: 55,
                             nSoilLevels: 4, nVertLevelsP1: 56)
Coordinates:
    cell                    (n_face) float64 6kB nan 1.0 2.0 ... 766.0 767.0
  * nSoilLevels             (nSoilLevels) float64 32B nan 1.0 2.0 3.0
  * nVertLevels             (nVertLevels) float64 440B nan 1.0 2.0 ... 53.0 54.0
  * nVertLevelsP1           (nVertLevelsP1) float64 448B nan 1.0 ... 54.0 55.0
  * time                    (time) datetime64[ns] 3kB 2020-01-20 ... 2020-03-01
Dimensions without coordinates: n_face
Data variables: (12/19)
    air_temperature         (time, n_face, nVertLevels) float32 56MB dask.array<chunksize=(1, 768, 55), meta=np.ndarray>
    h_oml                   (time, n_face) float32 1MB dask.array<chunksize=(1, 768), meta=np.ndarray>
    hu_oml                  (time, n_face) float32 1MB dask.array<chunksize=(1, 768), meta=np.ndarray>
    hv_oml                  (time, n_face) float32 1MB dask.array<chunksize=(1, 768), meta=np.ndarray>
    pressure                (time, n_face, nVertLevels) float32 56MB dask.array<chunksize=(1, 768, 55), meta=np.ndarray>
    relhum                  (time, n_face, nVertLevels) float32 56MB dask.array<chunksize=(1, 768, 55), meta=np.ndarray>
    ...                      ...
    tslb                    (time, n_face, nSoilLevels) float32 4MB dask.array<chunksize=(1, 768, 4), meta=np.ndarray>
    uReconstructMeridional  (time, n_face, nVertLevels) float32 56MB dask.array<chunksize=(1, 768, 55), meta=np.ndarray>
    uReconstructZonal       (time, n_face, nVertLevels) float32 56MB dask.array<chunksize=(1, 768, 55), meta=np.ndarray>
    vegfra                  (time, n_face) float32 1MB dask.array<chunksize=(1, 768), meta=np.ndarray>
    w                       (time, n_face, nVertLevelsP1) float32 57MB dask.array<chunksize=(1, 768, 56), meta=np.ndarray>
    xice                    (time, n_face) float32 1MB dask.array<chunksize=(1, 768), meta=np.ndarray>
ux_snow = uxds['snow']
ux_snow
<xarray.UxDataArray 'snow' (time: 329, n_face: 768)> Size: 1MB
dask.array<open_dataset-snow, shape=(329, 768), dtype=float32, chunksize=(1, 768), chunktype=numpy.ndarray>
Coordinates:
    cell     (n_face) float64 6kB nan 1.0 2.0 3.0 ... 764.0 765.0 766.0 767.0
  * time     (time) datetime64[ns] 3kB 2020-01-20 ... 2020-03-01
Dimensions without coordinates: n_face
Attributes:
    units:                      kg m^{-2}
    long_name:                  snow water equivalent
    healpix_nside:              8
    healpix_dim_name:           cell
    original_mpas_spatial_dim:  nCells
%%time
projection = ccrs.Robinson(central_longitude=-135.5808361)

ux_snow.mean('time').plot(
    projection=projection,
    cmap="Blues",
    features=["borders", "coastline"],
    title="Global snow: Jan-March mean",
    width=700,
)
CPU times: user 4.9 s, sys: 184 ms, total: 5.09 s
Wall time: 1.03 s
WARNING:param.GeoOverlayPlot02427: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.